home *** CD-ROM | disk | FTP | other *** search
/ El Mac 9 / El Mac 9.iso / Shareware / Applications / MathPad 2.4 / XFuns / XFun kit / linfit src / linfitPPC.c < prev   
Encoding:
C/C++ Source or Header  |  1996-02-13  |  1.6 KB  |  77 lines  |  [TEXT/CWIE]
  1. /*
  2.  linfit(XYpoints)
  3.  Fit data to a straight line using linear regression formula.
  4. */
  5.  
  6. /* This file implements linfit() as a native powerPC code fragment.
  7.    For comparison see:
  8.    linfit68K.c   -- 68K code seg without A4 globals
  9.    linfit68Kg.c  -- 68K code seg with A4 globals
  10. */
  11.  
  12. #include "callbackg.h"
  13. #include <fp.h>        // needed for sqrt()
  14.  
  15. #include <CodeFragments.h>
  16. /* Fragment initialization entry point. Must match symbol name in linker prefs */
  17. OSErr InitCFun(CFragInitBlockPtr iblk);
  18.  
  19. static short linfit(double *retval)
  20. {
  21.    EXPR arr;
  22.    double x,y,*iptr,*jptr,intercept;
  23.    double sumx,sumxx,sumy,sumyy,sumxy,Sxx,Sxy;
  24.    short isarray;
  25.    long n,count;
  26.    
  27.    sumx=0; sumxx=0; sumy=0; sumyy=0; sumxy=0;
  28.    n = 0;
  29.    MakeParmExpr(0,&arr);
  30.    ProbeExpr(arr,&x,&isarray,&count);
  31.    if(!isarray || count<2)
  32.    {
  33.     ErrMsg(" linfit() expects array of {x,y}",0L);
  34.     FreeExpr(arr);
  35.     return(FALSE);
  36.    }
  37.    AddIndex(&arr,&iptr);    // arr[i] is {x,y}
  38.    AddIndex(&arr,&jptr);    // j picks x or y
  39.    while(count--)
  40.    {
  41.     *jptr = 1;
  42.     if(EvalExpr(arr,&x) && x==x)
  43.     {
  44.      *jptr = 2;
  45.      if(EvalExpr(arr,&y) && y==y)
  46.      {
  47.       sumx += x;
  48.       sumxx += x*x;
  49.       sumy += y;
  50.       sumyy += y*y;
  51.       sumxy += x*y;
  52.       n++;
  53.      }
  54.     }
  55.     *iptr += 1;
  56.    }
  57.    Sxx=n*sumxx-sumx*sumx;
  58.    Sxy=n*sumxy-sumx*sumy;
  59.  
  60.    SetVarVal("slope",Sxy/Sxx);
  61.    intercept=(sumxx*sumy-sumx*sumxy)/(n*sumxx-sumx*sumx);
  62.    SetVarVal("intercept",intercept);
  63.    *retval = Sxy/sqrt(Sxx*(n*sumyy-sumy*sumy));        // return correlation coeff
  64.    
  65.    FreeExpr(arr);
  66.    return(TRUE);
  67. }
  68.  
  69.  
  70.  
  71. OSErr InitCFun(CFragInitBlockPtr iblk)
  72. {
  73.    AddCFun("linfit","XYpoints",&linfit,NULL);
  74.    return noErr;
  75. }
  76.  
  77.